
################################################################################
################################################################################
###                                                                          ###
### --- Filename: 2. Main Analysis.R --------------------------------------- ###
###                                                                          ###
### --- Description: This is the R script for the main analysis for -------- ###
### --- the paper "Economic Inequality and Public Support for -------------- ###
### --- Redistribution in Europe: A Cross-Sectional and Longitudinal ------- ###
### --- Multilevel Analysis" submitted to the IJPOR ------------------------ ###
###                                                                          ###
### --- Author: Valentin Velev --------------------------------------------- ###
###                                                                          ###
### --- Date: 11.12.2023 --------------------------------------------------- ###
###                                                                          ###
################################################################################
################################################################################

################################################################################
### --- 1. Preparation ----------------------------------------------------- ###
################################################################################

# install and load packages
required_packages <- c(
  # packages for loading and wrangling data
  "tidyverse", "haven", "readxl", "sjmisc", "sjlabelled", "fastDummies",
  # packages for (multiple) imputation
  "naniar", "imputeTS", "jomo", "mitml",
  # packages for statistical models
  "lme4", "lmerTest", "broom.mixed",
  # packages for post-estimation
  "influence.ME", "easystats", "sjPlot", "jtools", "interplot", "interactions",
  # other packages
  "ggpubr", "qqplotr", "NCmisc", "rstatix"
)

for (p in required_packages){
  
  if(!require(p, character.only=T)) install.packages(p, dependencies=T)
  library(p, character.only=T)
  
}

rm(p, required_packages)

# set seed to ensure replicability
set.seed(2946)

# global options
options(digits=6, scipen=999, max.print=10000)

# load necessary environment objects
load("Environment/reg_dat.RData")

################################################################################
### --- 2. Estimating mixed-effects models --------------------------------- ###
################################################################################

# create empty/null model (Model M0)
m0 <- lmer(redist ~
             1+(1|country)+(1|country:year), 
           data=reg_dat, REML=T)

# adding individual-level variables and year dummies (Model M1)
m1 <- lmer(redist ~
             income_cgm+lrscale_cgm+age_cgm+I(age_cgm^2)+male+edulvl+employment+
             married+badhealth+unionmember+year+(1|country)+(1|country:year),
           data=reg_dat, REML=T)

# adding BE&WE of Gini index (Model M2)
m2 <- lmer(redist ~
             income_cgm+lrscale_cgm+age_cgm+I(age_cgm^2)+male+edulvl+employment+
             married+badhealth+unionmember+gini_disp_be+gini_disp_we+year+
             (1|country)+(1|country:year),
           data=reg_dat, REML=T)

# adding BE&WE of GDP/capita (Model M3)
m3 <- lmer(redist ~
             income_cgm+lrscale_cgm+age_cgm+I(age_cgm^2)+male+edulvl+employment+
             married+badhealth+unionmember+gini_disp_be+gini_disp_we+gdp_pc_be+
             gdp_pc_we+year+(1|country)+(1|country:year),
           data=reg_dat, REML=T)

# adding random slopes for income (Model M4)
m4 <- lmer(redist ~
             income_cgm+lrscale_cgm+age_cgm+I(age_cgm^2)+male+edulvl+employment+
             married+badhealth+unionmember+gini_disp_be+gini_disp_we+gdp_pc_be+
             gdp_pc_we+year+(1+income_cgm|country)+(1+income_cgm|country:year),
           data=reg_dat, REML=T,
           control=lmerControl(optimizer="bobyqa",
                               optCtrl=list(maxfun=1e6),
                               calc.derivs=F))

# adding cross-level interactions (Gini BE&WE*income) (Model M5)
m5 <- lmer(redist ~
             income_cwc+lrscale_cgm+age_cgm+I(age_cgm^2)+male+edulvl+employment+
             married+badhealth+unionmember+gini_disp_be+gini_disp_we+gdp_pc_be+
             gdp_pc_we+(gini_disp_be*income_cwc)+(gini_disp_we*income_cwc)+
             year+(1+income_cwc|country)+(1+income_cwc|country:year),
           data=reg_dat, REML=T,
           control=lmerControl(optimizer="bobyqa",
                               optCtrl=list(maxfun=1e6),
                               calc.derivs=F))

# display results in regression table (Table 1 in Paper)
tab_model(m0, m1, m2, m3, m4, m5,
          title="Main Models - REML & Satterthwaite approx. (ESS 2002-2018)",
          dv.labels=c("M0", "M1", "M2", "M3", "M4", "M5"),
          string.est="b/SE", p.val="satterthwaite", p.style="stars", show.aic=T,
          show.se=T, show.ci=F, show.dev=T, collapse.se=T, digits=4, show.icc=T,
          digits.re=4, digits.rsq=4, file="Regression Tables/tab_1.doc")

tab_model(m0, m1, m2, m3, m4, m5,
          title="Main Models - REML & Satterthwaite approx. (ESS 2002-2018)",
          dv.labels=c("M0", "M1", "M2", "M3", "M4", "M5"),
          string.est="b/SE", p.val="satterthwaite", p.style="numeric", show.aic=T,
          show.se=T, show.ci=F, show.dev=T, collapse.se=T, digits=4, show.icc=T,
          digits.re=4, digits.rsq=4)

# display variance and covariance components for each model
data.frame(VarCorr(m0))
data.frame(VarCorr(m1))
data.frame(VarCorr(m2))
data.frame(VarCorr(m3))
data.frame(VarCorr(m4))
data.frame(VarCorr(m5))

# AIC & BIC values
jtools::summ(m0, digits=getOption("jtools-digits", default=2))
jtools::summ(m1, digits=getOption("jtools-digits", default=2))
jtools::summ(m2, digits=getOption("jtools-digits", default=2))
jtools::summ(m3, digits=getOption("jtools-digits", default=2))
jtools::summ(m4, digits=getOption("jtools-digits", default=2))
jtools::summ(m5, digits=getOption("jtools-digits", default=2))

################################################################################

# re-estimating Model M3 with standardized continuous variables
jtools::summ(m3, scale=T, n.sd=2, scale.only=T, digits=getOption("jtools-digits", default=4))

################################################################################
### --- 3. Interaction plots ----------------------------------------------- ###
################################################################################

# interaction plot 1: Gini BE * income CWC
cme_plot1 <- interplot(m5, var1="income_cwc", var2="gini_disp_be", hist=T, ralpha=.5, ci=.95) +
  geom_hline(yintercept=0, linetype="dashed", linewidth=.3) +
  labs(x="Gini Index [BE]", y="CME of Household Income") +
  theme(axis.title.x = element_text(size=9),
        axis.title.y = element_text(size=9),
        axis.text = element_text(size=7),
        plot.background = element_blank(),
        panel.background = element_rect(fill="white", color="black"),
        panel.grid.minor = element_blank(),
        panel.grid.major = element_blank(),
        text = element_text(family="serif", face="bold"),
        legend.background = element_rect(fill = "white", color = NA),  # White legend background
        legend.key = element_rect(fill = "white", color = NA)          # White legend key background
        )

# interaction plot 2: Gini WE * income CWC
cme_plot2 <- interplot(m5, var1="income_cwc", var2="gini_disp_we",hist=T, ralpha=.5, ci=.95) +
  geom_hline(yintercept=0, linetype="dashed", linewidth=.3) +
  labs(x="Gini Index [WE]", y="CME of Household Income") +
  theme(axis.title.x = element_text(size=9),
        axis.title.y = element_text(size=9),
        axis.text = element_text(size=7),
        plot.background = element_blank(),
        panel.background = element_rect(fill="white", color="black"),
        panel.grid.minor = element_blank(),
        panel.grid.major = element_blank(),
        text = element_text(family="serif", face="bold"),
        legend.background = element_rect(fill = "white", color = NA),  # White legend background
        legend.key = element_rect(fill = "white", color = NA)          # White legend key background
        )

# append plots (create Figure 1)
fig1 <- ggarrange(cme_plot1, cme_plot2, ncol=2)

# save plot as .png file
ggsave(filename="Figures/fig1.jpeg", plot=fig1, width=8, height=3.5, units="in", dpi=1750)

# save plot as PDF file
ggsave(filename="Figures/fig1.pdf", plot=fig1, width=8, height=3.5, units="in", dpi=1750, device=cairo_pdf)

# save plot as tiff file
tiff('Figures/fig1.tiff', width=8, height=3.5, units="in", res=1750, compression="lzw")
print(fig1)
dev.off()

################################################################################
### --- 4. Robustness checks ----------------------------------------------- ###
################################################################################

################################################################################
### --- 4.1 Jackknifing ---------------------------------------------------- ###
################################################################################

# create for loop objects
countries <- c("AT", "BE", "BG", "HR", "CY", "CZ", "DK", "EE", "FI", "FR", "DE",
               "GR", "HU", "IS", "IE", "IT", "LV", "LT", "NL", "NO", "PL", "PT",
               "RO", "SK", "SI", "ES", "SE", "CH", "UK")

reg_list = list()
counter = ""

coefficients_gini_be = data.frame()
standard_errors_gini_be = data.frame()
p_values_gini_be = data.frame()
coefficients_gini_we = data.frame()
standard_errors_gini_we = data.frame()
p_values_gini_we = data.frame()

# for loop
for (i in countries){
  
  counter = i
  
  # Model M3
  reg_list[[counter]] <- 
    summary(lmerTest::lmer(redist ~
                             income_cgm+lrscale_cgm+age_cgm+I(age_cgm^2)+male+
                             edulvl+employment+married+badhealth+unionmember+
                             gini_disp_be+gini_disp_we+gdp_pc_be+gdp_pc_we+
                             year+(1|country)+(1|country:year),
                           data=reg_dat[reg_dat$country!=i,], REML=T,
                           control=lmerControl(optimizer="bobyqa",
                                               optCtrl=list(maxfun=1e6),
                                               calc.derivs=F)))[["coefficients"]]
  
  # Gini BE estimates
  coefficients_gini_be <- 
    rbind(coefficients_gini_be, reg_list[[i]]["gini_disp_be","Estimate"])
  
  standard_errors_gini_be <- 
    rbind(standard_errors_gini_be, reg_list[[i]]["gini_disp_be","Std. Error"])
  
  p_values_gini_be <- 
    rbind(p_values_gini_be, reg_list[[i]]["gini_disp_be","Pr(>|t|)"])
  
  # Gini WE estimates
  coefficients_gini_we <- 
    rbind(coefficients_gini_we, reg_list[[i]]["gini_disp_we","Estimate"])
  
  standard_errors_gini_we <- 
    rbind(standard_errors_gini_we, reg_list[[i]]["gini_disp_we","Std. Error"])
  
  p_values_gini_we <- 
    rbind(p_values_gini_we, reg_list[[i]]["gini_disp_we","Pr(>|t|)"])
  
}

# create df with results of jackknife
jackknifed_vals <- cbind(coefficients_gini_be, standard_errors_gini_be, p_values_gini_be,
                         coefficients_gini_we, standard_errors_gini_we, p_values_gini_we) %>%
  setNames(c("coefficients_gini_be", "standard_errors_gini_be", "p_values_gini_be",
             "coefficients_gini_we", "standard_errors_gini_we", "p_values_gini_we")) %>%
  add_column(country=c("Austria", "Belgium", "Switzerland", "Croatia", "Cyprus", "Czechia",
                       "Denmark", "Estonia", "Finland", "France", "Germany", "Greece",
                       "Hungary", "Iceland", "Ireland", "Italy", "Latvia", "Lithuania",
                       "Netherlands", "Norway", "Poland", "Portugal", "Romania", "Slovakia",
                       "Slovenia", "Spain", "Sweden", "Switzerland", "United Kingdom"), 
             .before="coefficients_gini_be")

# remove unnecessary objects from environment
rm(coefficients_gini_be, standard_errors_gini_be, p_values_gini_be,
   coefficients_gini_we, standard_errors_gini_we, p_values_gini_we,
   countries, counter, i, reg_list)

# average values of estimates
# Gini BE
round(mean(jackknifed_vals$coefficients_gini_be),4)
round(mean(jackknifed_vals$standard_errors_gini_be),4)
round(mean(jackknifed_vals$p_values_gini_be),4)

# Gini WE
round(mean(jackknifed_vals$coefficients_gini_we),4)
round(mean(jackknifed_vals$standard_errors_gini_we),4)
round(mean(jackknifed_vals$p_values_gini_we),4)

################################################################################
### --- 4.2 LMMs with additional country-level variables ------------------- ###
################################################################################

# Model M3 with electoral system and BE&WE of vote share of social democrats
m3.1 <- lmer(redist ~
               income_cgm+lrscale_cgm+age_cgm+I(age_cgm^2)+male+edulvl+employment+
               married+badhealth+unionmember+
               gini_disp_be+gini_disp_we+gdp_pc_be+gdp_pc_we+
               elec_sys+left_gov_be+left_gov_we+
               year+(1|country)+(1|country:year),
             data=reg_dat, REML=T)

# Model M3 with BE&WE of unemployment variable
m3.2 <- lmer(redist ~
               income_cgm+lrscale_cgm+age_cgm+I(age_cgm^2)+male+edulvl+employment+
               married+badhealth+unionmember+
               gini_disp_be+gini_disp_we+gdp_pc_be+gdp_pc_we+
               unemp_be+unemp_we+
               year+(1|country)+(1|country:year),
             data=reg_dat, REML=T)

# Model M3 with BE&WE of immigration variable
m3.3 <- lmer(redist ~
               income_cgm+lrscale_cgm+age_cgm+I(age_cgm^2)+male+edulvl+employment+
               married+badhealth+unionmember+
               gini_disp_be+gini_disp_we+gdp_pc_be+gdp_pc_we+
               immig_be+immig_we+
               year+(1|country)+(1|country:year),
             data=reg_dat, REML=T)

# Model M3 with BE&WE of globalization index
m3.4 <- lmer(redist ~
               income_cgm+lrscale_cgm+age_cgm+I(age_cgm^2)+male+edulvl+employment+
               married+badhealth+unionmember+
               gini_disp_be+gini_disp_we+gdp_pc_be+gdp_pc_we+
               global_index_be+global_index_we+
               year+(1|country)+(1|country:year),
             data=reg_dat, REML=T)

# Model M3 with BE&WE of social expenditure
m3.5 <- lmer(redist ~
               income_cgm+lrscale_cgm+age_cgm+I(age_cgm^2)+male+edulvl+employment+
               married+badhealth+unionmember+
               gini_disp_be+gini_disp_we+gdp_pc_be+gdp_pc_we+
               soc_exp_be+soc_exp_we+
               year+(1|country)+(1|country:year),
             data=reg_dat, REML=T)

# Model M3 with all additional country-level variables
m3.6 <- lmer(redist ~
               income_cgm+lrscale_cgm+age_cgm+I(age_cgm^2)+male+edulvl+employment+
               married+badhealth+unionmember+
               gini_disp_be+gini_disp_we+gdp_pc_be+gdp_pc_we+
               elec_sys+left_gov_be+left_gov_we+unemp_be+unemp_we+immig_be+immig_we+
               global_index_be+global_index_we+soc_exp_be+soc_exp_we+
               year+(1|country)+(1|country:year),
             data=reg_dat, REML=T)

################################################################################

# create dataframe with coefficients from Model M3.6
coefs_m3.6 <- 
  broom.mixed::tidy(m3.6, conf.int=T, conf.level=.95)[,c("term", "estimate", "std.error",
                                                         "statistic", "conf.low", "conf.high")] %>%
  sjmisc::rename_variables(term=variable, statistic=tvalue, conf.low=ci95_low, conf.high=ci95_high)

# create coefficient plot of country-level variables from Model M3.6 (Figure 2)
fig2 <- ggplot(data=coefs_m3.6[c(16:17, 22:31),],
               aes(x=estimate, y=variable)) +
  geom_vline(xintercept=0, color="black", linetype="dashed", linewidth=.3) +
  geom_point(aes(x=estimate, y=variable), size=1.5, color="black") +
  geom_linerange(aes(y=variable, xmin=ci95_low, xmax=ci95_high), 
                 linewidth=.5, color="black") +
  scale_y_discrete(limits=c("soc_exp_we","soc_exp_be","global_index_we","global_index_be",
                            "immig_we","immig_be","unemp_we","unemp_be","left_gov_we",
                            "left_gov_be","gini_disp_we","gini_disp_be"),
                   labels=c("soc_exp_we"="Social Expenditure [WE]",
                            "soc_exp_be"="Social Expenditure [BE]",
                            "global_index_we"="Globalization Index [WE]",
                            "global_index_be"="Globalization Index [BE]",
                            "immig_we"="Prop. of Int. Migrants [WE]",
                            "immig_be"="Prop. of Int. Migrants [BE]",
                            "unemp_we"="Unemployment [WE]",
                            "unemp_be"="Unemployment [BE]",
                            "left_gov_we"="Vote Share of Soc. Dem. [WE]",
                            "left_gov_be"="Vote Share of Soc. Dem. [BE]",
                            "gini_disp_we"="Gini Index [WE]",
                            "gini_disp_be"="Gini Index [BE]")) +
  labs(x="Coefficient (\u03b2)", y="Variable") +
  theme(axis.title.x = element_text(size=8),
        axis.title.y = element_text(size=9),
        axis.text = element_text(size=6),
        plot.background = element_blank(),
        panel.background = element_rect(fill="white", colour="black"),
        panel.grid.minor = element_blank(),
        panel.grid.major = element_blank(),
        text = element_text(family="serif", face="bold"))

# save plot as .png file
ggsave(filename="Figures/fig2.png", plot=fig2, width=6, height=4, units="in", dpi=1500)

# save plot as PDF file
ggsave(filename="Figures/fig2.pdf", plot=fig2, width=6, height=4, units="in", dpi=1500, device=cairo_pdf)

# save plot as tiff file
tiff('Figures/fig2.tiff', width=6, height=3.5, units="in", res=1750, compression="lzw")
print(fig2)
dev.off()

################################################################################
### --- 4.3 LMMs with market inequality ------------------------------------ ###
################################################################################

# Model M3 with gini_mrkt_be and gini_mrkt_we
m3_alt <- lmer(redist ~
                 income_cgm+lrscale_cgm+age_cgm+I(age_cgm^2)+male+edulvl+employment+
                 married+badhealth+unionmember+
                 gini_mrkt_be+gini_mrkt_we+gdp_pc_be+gdp_pc_we+
                 year+(1|country)+(1|country:year),
               data=reg_dat, REML=T)

# Model M5 with gini_mrkt_be and gini_mrkt_we
m5_alt <- lmer(redist ~
                 income_cwc+lrscale_cgm+age_cgm+I(age_cgm^2)+male+edulvl+employment+
                 married+badhealth+unionmember+
                 gini_mrkt_be+gini_mrkt_we+gdp_pc_be+gdp_pc_we+
                 (gini_mrkt_be*income_cwc)+(gini_mrkt_we*income_cwc)+
                 year+(1+income_cwc|country)+(1+income_cwc|country:year),
               data=reg_dat, REML=T,
               control=lmerControl(optimizer="bobyqa",
                                   optCtrl=list(maxfun=1e6),
                                   calc.derivs=F))

# display results in regression table
tab_model(m3_alt, m5_alt,
          title="Model M3 and M5 - Market Gini index", dv.labels=c("M3", "M5"),
          string.est="b/SE", p.val="satterthwaite", p.style="numeric", show.aic=T,
          show.se=T, show.ci=F, show.dev=T, collapse.se=T, digits=4, show.icc=F,
          digits.re=4, digits.rsq=4, digits.p=4)

################################################################################
### --- 5. Save ------------------------------------------------------------ ###
################################################################################

save(m0, m1, m2, m3, m4, m5, file="Environment/main_models.RData")
save(m3.1, m3.2, m3.3, m3.4, m3.5, m3.6, file="Environment/additional_models_1.RData")
save(m3_alt, m5_alt, file="Environment/additional_models_2.RData")

